function res = GetXe131T1Rho(obj, duration, varargin)
    p=inputParser;
    p.addRequired('duration', @(x) x>=0);
    p.addParameter('amplitude', 1.0, @(x) x>0.0 && x<1.0);
    p.addParameter('detuning', 0.0, @isnumeric);
    p.addParameter('bandwidth', 5.0, @(x) x>0.0);
    p.addParameter('preDrivingTime', 10.0, @(x) x>0);
    p.addParameter('calibrationField', 0.02, @(x) x>0.0);
    p.addParameter('rate', 50, @(x) x>0);
    p.addParameter('isReport', false, @islogical);
    p.addParameter('description', 'normal', @ischar);
    p.parse(duration, varargin{:});
    
    try
        freq131 = obj.Xe131FIDSignal.fft.center;
    catch
        error('Cannot get Central frequency.');
    end
    
    tm = timer('ExecutionMode', 'singleShot', 'StartDelay', p.Results.preDrivingTime);
    tm.TimerFcn = @(myTimerObj, thisEvent) drvXe(obj, myTimerObj, freq131 + p.Results.detuning, p.Results.amplitude, p.Results.calibrationField);
    start(tm);
    
    
    rbSignal = obj.setRFMagnetometer('field', obj.MainField.bz0, 'bandwidth', p.Results.bandwidth);
    rbSignal.subscribe_signal(ziDemodTime.R, [2 1 1]);
    rbSignal.subscribe_signal(ziDemodTime.Theta, [2 1 2]);
    res.signal = rbSignal.getDemodSignal(p.Results.duration, 'rate', p.Results.rate);
    
    res.amplitude = p.Results.amplitude;
    res.detuning = p.Results.detuning;
    res.calibration = p.Results.calibrationField;
    res.timestamp = datetime('now');
    res.toc = tm.UserData;
    
    calibrationCoeff = obj.MainField.f0/4.7/obj.MainField.bz0*p.Results.calibrationField; %Rb85
    [res.calibrationValue, res.data_t1rho] = findTrigger(res, calibrationCoeff);
    [res.fitting, res.envelope] = fitenvelopt(res.data_t1rho, -1);
    res.XeField = res.fitting.result.a;
    res.T1Rho = res.fitting.result.tau;
    
    delete(tm);
    obj.by131.disableDriver();
    
    obj.Xe131T1Rho = res;
    fig = obj.PlotXe131T1Rho();
    %%
    if p.Results.isReport
        getReporter;
        rpt.NewRecord(['Xe131T1Rho_' p.Results.description]);
        rpt.AppendFigure(fig, sprintf('Xe131T1Rho %s', p.Results.description));
    end
    
    getSession;
    sess.addData('Xe131T1Rho', res, p.Results.description);
end

function drvXe(obj, timerObj, freq, amp, calibration)
    tic;
    obj.bz0.setDC( obj.MainField.bz0 + calibration); pause(2.0);
    t1 = toc;
    obj.bz0.setDC( obj.MainField.bz0); pause(1.0);
    t2 = toc;
    obj.bz0.setDC( obj.MainField.bz0 - calibration); pause(2.0);
    t3 = toc;
    obj.bz0.setDC( obj.MainField.bz0); pause(1.0);
    t4 = toc;
    obj.by131.setAC('frequency', freq, 'amplitude', amp, 'range', 10.0); 
    timerObj.UserData = [t1, t2, t3, t4];
end

function [calibrationValue, data_t1rho] = findTrigger(res, calibrationField)
    tList0 = res.signal.getNodeData(2).x;
    data0 = res.signal.getNodeData(2).y;
    [~, idx] = findpeaks(diff(data0), 'Npeak', 2, 'SortStr', 'descend');
    t0 = tList0(min(idx)+1);
    tList = tList0 - t0;
    
    idxCalibration = tList>0 & tList<res.toc(1);
    calibrationValue = mean(data0(idxCalibration));
    
    idxT1rho = tList >= res.toc(4);
    data_t1rho.x = tList(idxT1rho) - res.toc(4);
    data_t1rho.y = data0(idxT1rho)/calibrationValue*calibrationField;
end

function [fitting, env] = fitenvelopt(data, sign)
    if nargin == 1
        sign = 1.0;
    end
    tList = data.x;
    env=envelope(sign*data.y, 3000);

    [xData, yData] = prepareCurveData( tList, env );

    % Set up fittype and options.
    ft = fittype( 'a*exp(-x/tau)+c', 'independent', 'x', 'dependent', 'y' );
    opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
    opts.Display = 'Off';
    opts.Lower = [0 0 0];
    opts.StartPoint = [90 40 100];

    % Fit model to data.
    [fitting.result, fitting.gof] = fit( xData, yData, ft, opts );
end